Sesión 7

Curso: Análisis Estadístico de Datos en Salud

Percy Soto-Becerra, M.D., M.Sc(c)

DIS – IETSI, EsSalud

2023-06-16

  https://github.com/psotob91

Estimación y Pruebas de Hipótesis para Variables Respuesta Numéricas: Tres o más grupos

Agenda

  1. Estimación y Pruebas de Hipótesis para Variables Respuesta Numéricas: Tres o más grupos

  2. Pruebas de hipótesis para variables respuesta categóricas

Una motivación para comenzar

  • Se desea comparar 3 tratamientos (A, B y C) para evaluar su efecto en el IMC.

  • Una opción es comparar los promedios de IMC al final de la intervención.

  • Las 3 pruebas t de Student con sus valores p se muestran:

    • A vs B (p = 0.032)

    • B vs C (p = 0.062)

    • A vs C (p = 0.047)

Una motivación para comenzar (cont.)

  • ¿Con un nivel de significancia de 5%, podemos concluir que existe asociación entre el tratamiento y el IMC? ¿Entre qué tratamientos el IMC promedio es diferente?

  • Asumiremos que controlar la probabilidad de error tipo 1 es una cuestión fundamental en este estudio.

    • Hay un debate filósfico/estadístico de en qué circunstancias este control no es apropiado.

    • Corregir por la probabilidad de error tipo 1 no es la única opción para obtener inferencias válidas.

Probabilidad de error tipo 1

  • Probabilidad de error tipo 1 es también conocido como nivel de significancia.

  • \(1 - Pr(\beta)\) es conocida como Potencia o Poder estadístico.

¡Alerta! A más pruebas, más error

  • Si cada prueba t de Student tiene un nivel de significancia prefijado de 5% (\(\alpha = 0.05\)), la probabilidad de error tipo 1 global (de todas las pruebas en conjunto) no necesariamente será de 0.05.

  • Asumamos que las tres hipótesis son independientes entre sí, entonces tendríamos lo siguiente:

A más pruebas, más error (cont.)

  • Pruebas de hipótesis para comparar 3 o más grupos permiten mantener la probabilidad de error tipo 1 global por debajo de un nivel pre-fijado.

Imagen tomada de: Lee S; et al. What is the proper way to apply the multiple comparison test? Korean Joournal of Anesthesiology 2018; 71(5): 353-360. https://www.ncbi.nlm.nih.gov/pmc/articles/PMC6193594/#!po=73.0769

Soluciones al problema de la probabilidad de error tipo 1

  • No hay consenso entre expertos.

  • Las agencias reguladoras (p. ej., EMA, FDA, IETSI) optan por el ajuste de comparaciones múltiples.

    • Una opción es usar una secuencia ANOVA, luego pruebas de comparaciones múltiples.
    • Otras soluciones, más populares que la secuencia de ANOVA, las discutiremos al final.

Pruebas globables que controlan la probabilidad de error tipo 1

  • Podemos probar hipótesis acerca de parámetros de 3 o más grupos.

  • El enfoque similar que para dos poblaciones, pero las hipótesis nula y alterna son las siguientes:

H0: Todas las poblaciones tienen iguales medias|proporciones|medianas|tasas de la variable respuesta son iguales.

HA: Al menos dos poblaciones tienen diferentes medias|proporciones|medianas|tasas de la variable respuesta.

Análisis de Varianza de una vía

  • También llamada: one-way ANOVA

  • Permite probar la hipótesis de diferencia de medias entre dos o más poblaciones.

    • Los usamos comúnmente para probar diferencia de medias entre 3+.

  • Nombre del método hace referencia a que analiza la variabilidad (“varianza”):

    • Variabilidad de la variable numérica dentro de cada grupo y entre los grupos.
    • Si se cumplen ciertos supuestos, las varianzas permiten comparar medias.

Hipótesis estadísticas del ANOVA one-way

  • H0: Todas las medias poblacionales son iguales entre los grupos.

  • Ha: Al menos dos medias poblacionales son diferentes entre los grupos.

Lógica de ANOVA (1/4)

Mismas varianzas entre grupos (\(Var_{EG} = \triangle{X} = 10\)) entre escenarios, pero diferente varianza intragrupo (\(VAR_{IG}\)) entre escenarios)

Lógica de ANOVA (2/4)

Lógica de ANOVA (3/4)

Diferentes varianzas entre grupos (\(Var_{EG}\)) entre escenarios, pero misma varianza intragrupo (\(VAR_{IG} = \triangle{X} = 15\)) entre escenarios)

Lógica de ANOVA (4/4)

En resumen: ANOVA

Distribución y estadístico F

Particionando la variabilidad en ANOVA one-way

ANOVA one-way en R

library(rstatix)
datos %>% 
  anova_test(weight ~ treat, detailed = TRUE)
ANOVA Table (type II tests)

  Effect    SSn      SSd DFn DFd     F     p p<.05   ges
1  treat 49.021 7171.016   2 100 0.342 0.711       0.007
  • La prueba ANOVA one-way no genera resultados puntuales, debemos generarlos con otra función:
datos %>% 
  group_by(treat) %>% 
  get_summary_stats(weight, type = "mean_sd")
# A tibble: 3 × 5
  treat   variable     n  mean    sd
  <fct>   <fct>    <dbl> <dbl> <dbl>
1 Placebo weight      33  60.4  6.40
2 Dosis 1 weight      33  62.1 10.4 
3 Dosis 2 weight      37  61.2  8.13
  • También es bueno visualizar los patrones mediante gráfico de cajas:
library(ggpubr)
datos %>% 
  ggboxplot(x = "treat", y = "weight")

Supuestos

  • Aleatorización

  • Indenpendencia de observaciones

  • Normalidad

    • Supuesto de normalidad condicional: \(Y_i \sim Normal\), es decir \(Y|x \sim Normal\)
library(ggpubr)
datos %>% 
  ggqqplot("weight", facet.by = "treat")

+ Una mejor forma de evaluar la normalidad es usando los residuos.

  • Homogeneidad de varianzas (*)

    • Se podría evaluar si las rectas de los gráficos qqplot normal son paralelas o aproximadamente paralelas:
datos %>% 
  ggqqplot("weight", color = "treat")

  • Otra opción es evaluar los residuales. Esto lo veremos en regresión lineal.

  • Es preferible asumir que no hay homocedaticidad y realizar un ajuste robusto a heterogeneidad de varianzas.

datos %>% 
  anova_test(weight ~ treat, detailed = TRUE, white.adjust = TRUE)
ANOVA Table (type II tests)

  Effect DFn DFd     F    p p<.05
1  treat   2 100 0.329 0.72
  • No valores extremos.

(*) Este supuesto puede flexibilizarse:

  1. Usando algún método robusto de estimación de varianza.

  2. Usando remuestreo (Prueba de permutación o Bootstrapping)

datos %>% 
  anova_test(lh ~ treat, detailed = TRUE)
ANOVA Table (type II tests)

  Effect     SSn      SSd DFn DFd     F    p p<.05   ges
1  treat 470.911 21640.35   2 103 1.121 0.33       0.021
  • La prueba ANOVA one-way no genera resultados puntuales, debemos generarlos con otra función:
datos %>% 
  group_by(treat) %>% 
  get_summary_stats(lh, type = "mean_sd")
# A tibble: 3 × 5
  treat   variable     n  mean    sd
  <fct>   <fct>    <dbl> <dbl> <dbl>
1 Placebo lh          34 11.6  17.7 
2 Dosis 1 lh          34 13.7  16.6 
3 Dosis 2 lh          38  8.65  7.75
  • También es bueno visualizar los patrones mediante gráfico de cajas:
datos %>% 
  ggboxplot(x = "treat", y = "lh")

Supuestos

  • Aleatorización

  • Indenpendencia de observaciones

  • Normalidad

    • Supuesto de normalidad condicional: \(Y_i \sim Normal\), es decir \(Y|x \sim Normal\)
library(ggpubr)
datos %>% 
  ggqqplot("lh", facet.by = "treat")

+ Una mejor forma de evaluar la normalidad es usando los residuos.

  • Homogeneidad de varianzas (*)

    • Se podría evaluar si las rectas de los gráficos qqplot normal son paralelas o aproximadamente paralelas:
datos %>% 
  ggqqplot("lh", color = "treat")

  • Otra opción es evaluar los residuales. Esto lo veremos en regresión lineal.

  • A veces es preferible asumir que no hay homocedaticidad y realizar un ajuste robusto a heterogeneidad de varianzas.

datos %>% 
  anova_test(lh ~ treat, detailed = TRUE, white.adjust = TRUE)
ANOVA Table (type II tests)

  Effect DFn DFd     F     p p<.05
1  treat   2 103 1.502 0.227
  • No valores extremos.

(*) Este supuesto puede flexibilizarse:

  1. Usando algún método robusto de estimación de varianza.

  2. Usando remuestreo (Prueba de permutación o Bootstrapping)

Kruskal Wallis test

datos %>% 
  kruskal_test(weight ~ treat)
# A tibble: 1 × 6
  .y.        n statistic    df     p method        
* <chr>  <int>     <dbl> <int> <dbl> <chr>         
1 weight   106    0.0347     2 0.983 Kruskal-Wallis
  • La prueba Kruskal Wallis no genera resultados puntuales, debemos generarlos con otra función:
datos %>% 
  group_by(treat) %>% 
  get_summary_stats(weight, type = "median_iqr")
# A tibble: 3 × 5
  treat   variable     n median   iqr
  <fct>   <fct>    <dbl>  <dbl> <dbl>
1 Placebo weight      33     59   8.4
2 Dosis 1 weight      33     60   9  
3 Dosis 2 weight      37     61  12
  • También es bueno visualizar las distribuciones
library(ggpubr)
datos %>% 
  ggboxplot(y = "weight", x = "treat")

Supuestos

  • Aleatorización

  • Indenpendencia de observaciones

  • Misma distribución salvo por la mediana:

    • Esto significa que Kruskal Wallis también es perjudicado en cierto modo por la heterogeneidad de varianzas porque afecta el supuesto de igualdad de distribuciones

    • Podemos usar gráfico de densidades:

library(ggpubr)
datos %>% 
  ggdensity(x = "weight", 
            y = "..density..", 
            fill = "treat", 
            color = "treat", 
            add = "median")

+ Podemos usar gráfico de violin: Combina cajas y densidad
library(ggpubr)
datos %>% 
  ggviolin(x = "treat", 
            y = "weight", 
            color = "treat", 
            add = "boxplot")

Comparaciones múltiples

  • Hay varios enfoques de comparaciones múltiples.

  • Algunas son pre-hoc y otras post-hoc.

  • Veremos solo un par para ejemplificar:

  • Para ANOVA one-way clásico:
aov(weight ~ treat, data = datos) %>% 
  tukey_hsd()
# A tibble: 3 × 9
  term  group1  group2 null.value estimate conf.low conf.high p.adj p.adj.signif
* <chr> <chr>   <chr>       <dbl>    <dbl>    <dbl>     <dbl> <dbl> <chr>       
1 treat Placebo Dosis…          0    1.72     -3.24      6.68 0.689 ns          
2 treat Placebo Dosis…          0    0.745    -4.08      5.57 0.928 ns          
3 treat Dosis 1 Dosis…          0   -0.974    -5.80      3.85 0.881 ns
  • Si homogeneidad de varianzas no se cumple:
datos %>% 
  games_howell_test(weight ~ treat)
# A tibble: 3 × 8
  .y.    group1  group2  estimate conf.low conf.high p.adj p.adj.signif
* <chr>  <chr>   <chr>      <dbl>    <dbl>     <dbl> <dbl> <chr>       
1 weight Placebo Dosis 1    1.72     -3.42      6.85 0.701 ns          
2 weight Placebo Dosis 2    0.745    -3.43      4.91 0.904 ns          
3 weight Dosis 1 Dosis 2   -0.974    -6.39      4.44 0.903 ns
  • En epidemiología e investigación clínica se prefiere ajustes más conservadores. Uno de ellos es el ajuste de Bonferroni (uno de los más conservadores)
library(emmeans)
datos %>% 
  emmeans_test(weight ~ treat, p.adjust.method = "bonferroni")
# A tibble: 3 × 9
  term  .y.    group1  group2     df statistic     p p.adj p.adj.signif
* <chr> <chr>  <chr>   <chr>   <dbl>     <dbl> <dbl> <dbl> <chr>       
1 treat weight Placebo Dosis 1   100    -0.824 0.412     1 ns          
2 treat weight Placebo Dosis 2   100    -0.367 0.714     1 ns          
3 treat weight Dosis 1 Dosis 2   100     0.480 0.632     1 ns          
datos %>% 
  dunn_test(weight ~ treat)
# A tibble: 3 × 9
  .y.    group1  group2     n1    n2 statistic     p p.adj p.adj.signif
* <chr>  <chr>   <chr>   <int> <int>     <dbl> <dbl> <dbl> <chr>       
1 weight Placebo Dosis 1    33    33    0.173  0.863     1 ns          
2 weight Placebo Dosis 2    33    37    0.148  0.882     1 ns          
3 weight Dosis 1 Dosis 2    33    37   -0.0300 0.976     1 ns

Pruebas de hipótesis para variables respuesta categóricas

Agenda

  1. Estimación y Pruebas de Hipótesis para Variables Respuesta Numéricas: Tres o más grupos

  2. Pruebas de hipótesis para variables respuesta categóricas

Tablas de contingencia

  • Las tablas de contingencia nos permiten describir la distribución de frecuencias de dos variables categóricas.

    • Proporciones marginales, fila y columna.
CrossTable(datos_fum$fumar, datos_fum$cancer)

 
   Cell Contents
|-------------------------|
|                       N |
| Chi-square contribution |
|           N / Row Total |
|           N / Col Total |
|         N / Table Total |
|-------------------------|

 
Total Observations in Table:  550 

 
                | datos_fum$cancer 
datos_fum$fumar |    Cancer | No cancer | Row Total | 
----------------|-----------|-----------|-----------|
           Fuma |        20 |        80 |       100 | 
                |    10.000 |     1.111 |           | 
                |     0.200 |     0.800 |     0.182 | 
                |     0.364 |     0.162 |           | 
                |     0.036 |     0.145 |           | 
----------------|-----------|-----------|-----------|
        No fuma |        35 |       415 |       450 | 
                |     2.222 |     0.247 |           | 
                |     0.078 |     0.922 |     0.818 | 
                |     0.636 |     0.838 |           | 
                |     0.064 |     0.755 |           | 
----------------|-----------|-----------|-----------|
   Column Total |        55 |       495 |       550 | 
                |     0.100 |     0.900 |           | 
----------------|-----------|-----------|-----------|

Prueba de hipótesis para una proporción

  • A menudo, queremos probar hipótesis u obtener estimados acerca de la proporción del valor de una variable categórica.

  • Por ejemplo, queremos probar la hipótesis de que la proporción de fumadores es mayor a 20% en la población de estudio.

datos_fum %>% 
  tabyl(fumar) 
   fumar   n   percent
    Fuma 100 0.1818182
 No fuma 450 0.8181818
  • Tenemos varias opciones para probar hipótesis o estimar intervalos de confianza:
  • Dos colas:
prop_test(x = 100, n = 550, p = 0.2, detailed = TRUE) %>% 
  gt()
n n1 estimate statistic p df conf.low conf.high method alternative p.signif
550 100 0.1818182 1.025568 0.311 1 0.1509869 0.2171871 Prop test two.sided ns
  • Cola superior:
prop_test(x = 100, n = 550, p = 0.2, alternative = "greater", 
          detailed = TRUE) %>% 
  gt()
n n1 estimate statistic p df conf.low conf.high method alternative p.signif
550 100 0.1818182 1.025568 0.844 1 0.1554933 1 Prop test greater ns
  • Dos colas:
binom_test(x = 100, n = 550, p = 0.2, detailed = TRUE) %>% 
  gt()
n estimate statistic p parameter conf.low conf.high method alternative p.signif
550 0.1818182 100 0.3111145 550 0.1504567 0.2166426 Exact binomial test two.sided ns
  • Cola superior:
binom_test(x = 100, n = 550, p = 0.2, detailed = TRUE) %>% 
  gt()
n estimate statistic p parameter conf.low conf.high method alternative p.signif
550 0.1818182 100 0.3111145 550 0.1504567 0.2166426 Exact binomial test two.sided ns

Pruebas de hipótesis para dos variables categóricas

  • Otras veces podemos estar interesados en probar alguna hipótesis sobre la relación de dos variables categóricas.
datos_fum2 %>%
  tabyl(grupo, fumar)
   grupo Fuma No fuma
 Grupo 1   50      56
 Grupo 2  100      13
 Grupo 3  139      17
 Grupo 4   80      22
datos_fum2 %>%
  tabyl(grupo, fumar) %>% 
  select(-1) %>%
  row_wise_prop_test(detailed = TRUE) %>% 
  gt()
group n n1 n2 estimate1 estimate2 statistic p df conf.low conf.high method alternative p.adj p.adj.signif
1 477 50 56 0.1355014 0.5185185 68.71417683 1.14e-16 1 -0.48949836 -0.2765360 Prop test two.sided 4.56e-16 ****
2 477 100 13 0.2710027 0.1203704 9.66967143 1.87e-03 1 0.06834047 0.2329242 Prop test two.sided 3.74e-03 **
3 477 139 17 0.3766938 0.1574074 17.27140110 3.24e-05 1 0.12867391 0.3098988 Prop test two.sided 9.72e-05 ****
4 477 80 22 0.2168022 0.2037037 0.02515029 8.74e-01 1 -0.07970374 0.1059007 Prop test two.sided 8.74e-01 ns
  • Se crea la tabla de contingencia:
datos_fum2 %>% 
  tabyl(grupo, fumar) 
   grupo Fuma No fuma
 Grupo 1   50      56
 Grupo 2  100      13
 Grupo 3  139      17
 Grupo 4   80      22
  • Se elimina la primera columna desde la izquierda:
datos_fum2 %>% 
  tabyl(grupo, fumar) %>%
  select(-1)
 Fuma No fuma
   50      56
  100      13
  139      17
   80      22
  • Se aplica la prueba de hipótesis:
datos_fum2 %>% 
  tabyl(grupo, fumar) %>%
  select(-1) %>%
  chisq_test() %>%
  gt()
n statistic p df method p.signif
477 75.50793 2.82e-16 3 Chi-square test ****
  • Se pueden visualizar algunas medidas de interés como las frecuencias esperadas:
datos_fum2 %>% 
  tabyl(grupo, fumar) %>%
  select(-1) %>%
  chisq_test() %>%  
  chisq_descriptives() %>%
  gt()
Var1 Var2 observed prop row.prop col.prop expected resid std.resid
A Fuma 50 0.10482180 0.4716981 0.1355014 82.00000 -3.5338088 -8.4209793
B Fuma 100 0.20964361 0.8849558 0.2710027 87.41509 1.3460362 3.2382663
C Fuma 139 0.29140461 0.8910256 0.3766938 120.67925 1.6677351 4.2724900
D Fuma 80 0.16771488 0.7843137 0.2168022 78.90566 0.1231964 0.2920041
A No fuma 56 0.11740042 0.5283019 0.5185185 24.00000 6.5319726 8.4209793
B No fuma 13 0.02725367 0.1150442 0.1203704 25.58491 -2.4880440 -3.2382663
C No fuma 17 0.03563941 0.1089744 0.1574074 35.32075 -3.0826795 -4.2724900
D No fuma 22 0.04612159 0.2156863 0.2037037 23.09434 -0.2277190 -0.2920041
  • TAmbién se peude obtener con gmodels
CrossTable(datos_fum2$grupo, datos_fum2$fumar, chisq = TRUE)

 
   Cell Contents
|-------------------------|
|                       N |
| Chi-square contribution |
|           N / Row Total |
|           N / Col Total |
|         N / Table Total |
|-------------------------|

 
Total Observations in Table:  477 

 
                 | datos_fum2$fumar 
datos_fum2$grupo |      Fuma |   No fuma | Row Total | 
-----------------|-----------|-----------|-----------|
         Grupo 1 |        50 |        56 |       106 | 
                 |    12.488 |    42.667 |           | 
                 |     0.472 |     0.528 |     0.222 | 
                 |     0.136 |     0.519 |           | 
                 |     0.105 |     0.117 |           | 
-----------------|-----------|-----------|-----------|
         Grupo 2 |       100 |        13 |       113 | 
                 |     1.812 |     6.190 |           | 
                 |     0.885 |     0.115 |     0.237 | 
                 |     0.271 |     0.120 |           | 
                 |     0.210 |     0.027 |           | 
-----------------|-----------|-----------|-----------|
         Grupo 3 |       139 |        17 |       156 | 
                 |     2.781 |     9.503 |           | 
                 |     0.891 |     0.109 |     0.327 | 
                 |     0.377 |     0.157 |           | 
                 |     0.291 |     0.036 |           | 
-----------------|-----------|-----------|-----------|
         Grupo 4 |        80 |        22 |       102 | 
                 |     0.015 |     0.052 |           | 
                 |     0.784 |     0.216 |     0.214 | 
                 |     0.217 |     0.204 |           | 
                 |     0.168 |     0.046 |           | 
-----------------|-----------|-----------|-----------|
    Column Total |       369 |       108 |       477 | 
                 |     0.774 |     0.226 |           | 
-----------------|-----------|-----------|-----------|

 
Statistics for All Table Factors


Pearson's Chi-squared test 
------------------------------------------------------------
Chi^2 =  75.50793     d.f. =  3     p =  2.819929e-16

  • Cuando el número de celdas con frecuencais esperadas < 5 es mayor a 20% del total de celdas, la prueba Chi-2 no es válida.

  • En estos casos, se puede optar por la prueba Fisher como una alternativa más robusta.

  • Se crea la tabla de contingencia:

datos_fum2 %>% 
  tabyl(grupo, fumar) 
   grupo Fuma No fuma
 Grupo 1   50      56
 Grupo 2  100      13
 Grupo 3  139      17
 Grupo 4   80      22
  • Se elimina la primera columna desde la izquierda:
datos_fum2 %>% 
  tabyl(grupo, fumar) %>%
  select(-1)
 Fuma No fuma
   50      56
  100      13
  139      17
   80      22
  • Se aplica la prueba de hipótesis:
datos_fum2 %>% 
  tabyl(grupo, fumar) %>%
  select(-1) %>%
  fisher_test() %>%
  gt()
n p p.signif
477 6.1e-15 ****
  • También podemos usar gmodels:
CrossTable(datos_fum2$grupo, datos_fum2$fumar, fisher = TRUE)

 
   Cell Contents
|-------------------------|
|                       N |
| Chi-square contribution |
|           N / Row Total |
|           N / Col Total |
|         N / Table Total |
|-------------------------|

 
Total Observations in Table:  477 

 
                 | datos_fum2$fumar 
datos_fum2$grupo |      Fuma |   No fuma | Row Total | 
-----------------|-----------|-----------|-----------|
         Grupo 1 |        50 |        56 |       106 | 
                 |    12.488 |    42.667 |           | 
                 |     0.472 |     0.528 |     0.222 | 
                 |     0.136 |     0.519 |           | 
                 |     0.105 |     0.117 |           | 
-----------------|-----------|-----------|-----------|
         Grupo 2 |       100 |        13 |       113 | 
                 |     1.812 |     6.190 |           | 
                 |     0.885 |     0.115 |     0.237 | 
                 |     0.271 |     0.120 |           | 
                 |     0.210 |     0.027 |           | 
-----------------|-----------|-----------|-----------|
         Grupo 3 |       139 |        17 |       156 | 
                 |     2.781 |     9.503 |           | 
                 |     0.891 |     0.109 |     0.327 | 
                 |     0.377 |     0.157 |           | 
                 |     0.291 |     0.036 |           | 
-----------------|-----------|-----------|-----------|
         Grupo 4 |        80 |        22 |       102 | 
                 |     0.015 |     0.052 |           | 
                 |     0.784 |     0.216 |     0.214 | 
                 |     0.217 |     0.204 |           | 
                 |     0.168 |     0.046 |           | 
-----------------|-----------|-----------|-----------|
    Column Total |       369 |       108 |       477 | 
                 |     0.774 |     0.226 |           | 
-----------------|-----------|-----------|-----------|

 
Fisher's Exact Test for Count Data
------------------------------------------------------------
Alternative hypothesis: two.sided
p =  6.10405e-15

Hagamos una pausa

Tomemos un descanso de 5 minutos…

Estire las piernas …

Deje de ver las pantallas …

… cualquier , las del celular también.

05:00

¡Gracias!
¿Preguntas?




@psotob91

https://github.com/psotob91

percys1991@gmail.com

Análisis Estadístico de Datos en Salud 2023 - Día 7